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ABSTRACT 



Application of classic triangulation methods will allow 
the location of a radar to be determined by passive sensors. 
Through the use of modern digital signal processing techniques 
this estimate can be made in a simpler fashion using a 
conventional receiver. 

In this thesis a technique is developed for time 
difference of arrival (TDOA) estimation using a frequency 
domain based correlation detector driven by an envelope 
detector. Time lag boundaries are defined on the output of the 
correlator. A fixed detection threshold is calculated to 
permit constant false alarm rate (CFAR) detection. The 
performance of the correlation detector is plotted as a 
receiver operating characteristic (ROC) curve as a function of 
signal to noise ratio (SNR) . An interactive MATLAB software 
program is provided to perform either spectral domain or time 
domain based correlation. 

Spectral domain based correlation uses the Fast Fourier 
Transform (FFT) . Implicit with the use of the FFT are finite 
arithmetic internal processing errors which are modeled as 
independent uncorrelated noise sources. A method is presented 
to account for SNR degradation at the output of the FFT. 
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I . INTRODUCTION 



A . BACKGROUND 

Interception of radar transmissions by passive sensors is 
one of the responsibilities of electronic intelligence. 
Analysis of the signal of a radar can provide an estimate of 
the distance and direction of an emitter relative to an 
intercepting receiver. Knowing the location is an important 
consideration in the evaluation of strategic and tactical 
capabilities of the radar. 

Triangulation is a traditional technique where the angle 
of arrival of a radar signal at two spatially separated 
receivers provides an estimate of the transmitters location. 
With the advent of the computer, digital signal processing 
algorithms can be implemented that provide emitter range 
information and hence allow localization in a sequential 
sense. One of these processing algorithms is called the time 
difference of arrival (TDOA) method. 

B. TIME DIFFERENCE OF ARRIVAL 

The TDOA algorithm relies on two bistatic passive sensors 
attempting to receive the transmission of an active radar. 
Reception by the two sensors is preferentially done in the 
main lobe of the transmitting radar, but due to the 
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nonsynchronous nature of the task^ reception through one of 
the sidelobes of the transmitting radar is expected with a 
resulting decrease in received signal power. 

Given a transmitted radar pulse f(t)^ two physically 
separated passive sensors will most likely receive this pulse 
at times ti and ti+t. The time difference t is proportional to 
their differential distance relative to the emitter. This time 
difference of arrival can be used in the estimation of the 
differential range of the emitter to the receivers. The 
reception of the first pulse in a pulse burst by the sensors 
is paramount in performing the TDOA estimate. 

The concept of measuring the TDOA of a radar pulse between 
two separated receivers can be modeled by a hyperbola. The two 
sensors whose positions are known are located at the two focal 
points of the hyperbola. A transmitting radar is located in 
the area surrounding the two sensors. Figure 1 shows this 
model . 

A hyperbola is the locus of points in which the distance 
from any one point on the hyperbola to one focus differs by a 
constant amount from the distance to the other focus. This 
differential distance is proportional to a difference in 
arrival time X, of a radar pulse as detected between the two 
sensors. Once X is measured the hyperbola can be drawn. If the 
sensors (moving aircraft or nongeostationary satellites) 
repeat the differential arrival time measurements at several 
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time intervals, a series of hyperbolas will be generated, 
fixing the location of the transmitting radar. Alternately, 
the angle of arrival at either sensor can be used to define 
the position of the transmitter on the hyperbola. A good 
introduction into the use of the hyperbola in navigation is 
found in [Ref. l;p. 27]. 




Figure 1. Time difference of arrival model. 

This discussion asstimes that the radar signal received by 
the two sensors is induced by the direct wave from the 
transmitter. The effects of refraction and reflection on 
propagating radar waves are not covered here. 
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The purpose of this thesis is to develop and test an 
algorithm to estimate the differential arrival time x, of a 
pulsed radar signal collected by two passive sensors. If the 
pulse burst is collected in the time domain, time correlation 
can be used to generate a TDOA estimate. The Fast Fourier 
Transform (FFT) can also be used to perform the correlation of 
the two received signals (i.e., fast correlation). 

Additive noise n(t) , superimposed onto a transmitted radar 
signal f (t) by the environment and the radar receiver, creates 
a composite signal s (t) at the receiver 

s(t) = fit) ^ nit) . (1) 

The use of the Fourier Transform for TDOA estimation is 
attractive because of the processing gain (PG) a radar signal 
receives in being transformed from the time domain to the 
frequency domain during the detection process. FFT processing 
gain for real valued signals can be approximated by 

PG idB) = (log 2 [number points in data record] ~1)’3 . (2) 

See for example [Ref. 2:p. 33]. The signal to noise ratio 
(SNR) at the output of a FFT (SNRour) is a function of the PG 
and the SNR at the input to the FFT (SNRjn) • This relation is 
given by 
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(3) 



SNRq^j. - SNRj^i^ + PG . 

The corruptive effects of additive noise are reduced by the 
transformation, allowing the frequency domain detection of the 
signal . 

C . RADAR TYPES 

The TDOA algorithm is primarily applied to pulsed radars. 
The continuous wave (CW) radar provides target radial velocity 
through the Doppler shift of its return. The CW radar does not 
use a pulsed transmission. Therefore, the TDOA algorithm 
cannot not be directly applied to the reception of CW signals. 

Pulsed radars can be divided into two broad categories, 
ordinary low pulse repetition frequency (PRF) pulsed radars 
and the pulse Doppler radars. 

Low PRF radars can yield unambiguous range information, 
while their radar returns do not give target velocity 
information directly. Generally because of the long duty cycle 
of the low PRF radar, the probability of receiving all radar 
pulses in a pulse burst by two passive sensors is great. The 
TDOA algorithm can then estimate the differential distance to 
the emitter using the low PRF burst. 

The pulse Doppler radar uses coherent transmission and 
reception with a moderately high PRF [Ref. 3:p. 17.1]. It 
gives both range and radial velocity information about a 
target, though not with the same accuracy of the ordinary 
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pulsed radar or CW radar respectively. In the pulse Doppler 
radar, as in all non-bistatic radars, the receiver must be 
turned off while the transmitter is on. Because of this, a 
range ambiguity exists. Echoes having delay times equal to an 
integral multiple of the PRF will be undetected. Target echoes 
that remain within this blanked time interval will remain 
undetected. A second source of range ambiguity exists. For a 
target located at a distance r^ from the tracking radar, all 
other targets in the main beam of the radar a distance 



z^, 2i^, 3r^, 4r^, . . . , kz^ 
wheze k = integez . 



( 4 ) 



will appear to be approximately the same distance from the 
radar . 

Pulse Doppler radars can be divided into a low, medium and 
high PRF class. Both the low and medium PRF pulse Doppler 
radar transmissions can be passively intercepted and a 
relatively simple TDOA estimate can be formed. 

Reception of each pulse in a pulse burst from a high PRF 
pulse Doppler radar by two independent passive sensors can be 
difficult. If either sensor misses the first pulse, the TDOA 
estimate will be in error. This is true for any radar, but has 
a higher likelihood of occurring in the high PRF radars due to 
the shorter pulse widths used. Should the pulse burst be coded 
so that a nonuniform time pulse sequence is generated (i.e.. 
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staggered PRF) , then missed pulses by either sensor will not 
seriously degrade the validity of the TDOA estimate. 

D. FAST FOORIER TRANSFORM OUTPUT SIGNAL TO NOISE RATIO 

The Fourier transform is a computationally intensive 
mathematical operation. When implemented on an FFT processor 
using fixed point arithmetic, internal processing errors are 
generated. These errors are a function of the transform size, 
and the number of bits used internally in the FFT processor. 

The primary building block of the FFT is the butterfly 
algorithm. There are three dominant processing errors that 
occur within the butterfly. They are scaling, truncation and 
trigonometric errors. These can be modeled as independent, 
uncorrelated noise sources. The cumulative effect of these 
noise sources is to reduce the potential SNR at the output of 
the butterfly, and therefore, the FFT processor. 

This thesis will also examine the effects of noise caused 
by FFT processing, and will present a method to account for 
the degradation of the output SNR due to finite arithmetic. 
Certain FFT implementations will keep the signal energy larger 
with respect to the three noise power sources than others. 
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II . TIME DOMAIN BASED CORRELATION 



A. BACKGROUND 

The correlation function R(T) measures the degree of 
similarity between two signals s (t) and g(t) and is described 
by 

i?(x) = t) 5 t(t + t) dt 
where g{x + t) = g{C) displaced by the shift x . 

Correlation can be performed in a radar receiver between 
a target return corrupted by additive noise, and a replica of 
the transmitted signal maintained by the radar. This replica 
is designed into the frequency response of the matched filter 
of the radar receiver. [Ref. 4:p. 373]. 

Correlating two signals produces an output that does not 
resemble either of the two input signals. The shape of the 
received waveform is destroyed during the correlation process. 
The useful item of correlation is an amplitude peak, where the 
location along the time axis of this peak indicates the amount 
of time that one signal lags the other. 

If the correlation process occurs between a signal and a 
delayed replica of itself, the process is called 
autocorrelation. If two dissimilar signals are correlated, the 
process is called crosscorrelation. 
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B . AUTOCORRELATION 

Autocorrelation is the process by which two identical 
signals, one a delayed replica of the other, are correlated. 
The autocorrelation function of a discrete time sequence x(i) 
is defined by [Ref. 5:p. 556] 



w-i-m 

i=0 

where 1 = shift operator 

N = number of data points in x(n) 



( 6 ) 



This function is not scaled with respect to the shift operator 
nor the number of data points. Figure 2(a) shows the effect of 
autocorrelation on a burst of 13 pulses of unit amplitude. 
Normalized autocorrelation produces a triangular envelope 
waveform of height one at a time lag of zero. The location of 
this peak indicates that the time shift for maximum overlap 
between the pulse burst and a shifted replica of itself is 
zero . 

For illustration, a Barker coded sequence of pulses of 
length 13 is also autocorrelated in Figure 2 (b) . The Barker 
code is given by +++++ — ++-+-+, where + and - could represent 
0 and K radians carrier phase shift, respectively. Barker 
coding of a pulse sequence constructively modifies the 
correlation output to magnify the time lag peak. The zero lag 
point has the greatest amplitude. Adjacent peaks are 
significantly reduced in amplitude. Because of this reduction. 
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Figure 2. (a) Normalized autocorrelation of 13 equal 
amplitude pulses, (b) Normalized autocorrelation of 13 Barker 
coded pulses. 



if noise is added to the received signal, a Barker coded radar 
pulse burst will have a correlation peak that is not as easily 
confused with adjacent peaks. 

Now consider the autocorrelation of a random process. An 
independent, identically distributed zero mean Gaussian 
sequence with variance of one is shown in Figure 3. The 
autocorrelation of this sequence produces a peak at zero time 
lag. There is no strong correlation of the Gaussian noise 
except at zero time lag. The autocorrelation of Rayleigh 
distributed noise is shown in Figure 4. 
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Figure 4. (a) Rayleigh distributed noise sequence, 

(b) Normalized autocorrelation of Rayleigh noise. 
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Rayleigh distributed noise has a mean (dc) value of 



mean value - o\ — 
\ 2 



( 7 ) 



where o = standard deviation 



The autocorrelation function of Rayleigh noise has two 
components. A triangular envelope caused by the dc component 
of the Rayleigh noise, and a weighted delta function at the 
apex of the triangle which indicates the time shift of maximum 
overlap of the two sequences. Clearly maximum overlap occurs 
at a time shift of zero for any autocorrelation function. 

C. CROSSCORRELATION 

Crosscorrelation is the process wherein two different 
signals are correlated. The crosscorrelation function of a 
discrete time sequence is defined by 



This function is not scaled with respect to the shift operator 
nor the number of data points. Figure 5 shows the 
crosscorrelation of two uncorrelated Gaussian noise sequences 
and two uncorrelated Rayleigh noise sequences. Clearly, no 
dominant amplitude peak is visible in either crosscorrelation 
plot. Lack of a dominant peak indicates there is no strong 
correlation between the two sets of noise sequences. 







( 8 ) 



where x(i) and y(i) are two different sequences. 
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Figure 5. (a) Normalized Gaussian noise crosscorrelation, 

(b) Normalized Rayleigh noise crosscorrelation. 



D . CORREULTION COEFFICIENT 

A measure of the linear dependence between any two zero 
mean stationary time functions x(t) and y(t) is the 
correlation coefficient function and is defined [Ref. 6:p. 48] 
as 

(T) - 

where |p^(t) | i 1 for all t. 

This f\anction is the crosscorrelation fvinction normalized by 
the square root of the maximum values of the individual 
autocorrelation functions. The normalizing factor is not a lag 
dependent quantity. 
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Given two time series x(t) and y(t) where a is any 
positive number 

p^(-r) = 1 if x{t) = a-y(t) 

^ if x(t) = -a-yic) (10) 

P^(t) = 0 if x(t) and y(t) are uncorrelated. 



Another way of defining a normalized correlation 
coefficient for a discrete time series of finite duration is 
given by 



Pxx(-Z) 



W-l-UI 

x(i) yii + i) 

i=0 



N-1 



> i=0 






w-i-UI 



y(i+l)^ 



1=0 



( 11 ) 



E. SQUARE-LAW DETECTION AND CORRELATION 

If x(t) is a real valued function of time, the Hilbert 
transform of x(t) is defined by [Ref. 6:p. 484] as 



x( t) = H [x( t) ] 




x(u) 

X ( t - u) 



du 



= x(t) ♦ (^) 

X t 

where * = convolution operation . 



( 12 ) 



The output of a square-law envelope detector u(t), can be 
described by 
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( 13 ) 



u ( t) = (t) + it) 

where x{t) = real valued time function 

input to square-law detector 



The correlation coefficient function of two time functions 
X (t) and y(t) is defined in Equation 9. The Hilbert transform 
of the correlation function is defined as 



= H [p;g,(x) ] 

< 14 ) 

where and Oy = standard deviation of x and y 

respectively . 



The correlation coefficient for two square-law envelope 
detected signals, Pu^(T), is defined by [Ref. 6:p. 512] as 

Pu.(^> = P^(^> - P^(^) 



where u(t) = square-law detected signal of x{t) 
v{t) = square-law detected signal of y(t) 



( 15 ) 



The quantity p„(t) produces a sharper correlation peak than 
Pj 5 y(t) . This sharpening of the correlation function peak output 
can aid in locating the TDOA correlation peak. An example of 
both Pxy(x) and Puv('^) is plotted in Figure 6. 



F. ENVELOPE CORRELATION AND RELATED STATISTICS 

In this thesis, time domain based correlation will be 
performed on the output of an envelope detector. An envelope 
detector can be easily implemented in hardware using a diode. 



15 




Figure 6. (a) Normalized correlation coefficient, (b) Hilbert 
transform of the normalized correlation coefficient (dotted 
line) . 



To derive the expected value and variance of the 
crosscorrelation function, two cases must be considered. These 
cases are noise only present, and signal plus noise present in 
both channels of the correlator. Under the noise only case, 
the statistics for the correlator output are easily derived. 

We consider two real, independent identically distributed, 
zero mean noise series at the input to the correlator. Each 
series has zero mean and variance . It is shown in Appendix 
A, that the expected value of the output of the 
crosscorrelation function r,y(P) is 

E[r^yil) ] = 0 . ( 16 ) 
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The variance of the output of the crosscorrelation function 
r*y(8) is 



Since equal variances are assumed for both time series, is 
the square of the variance of either noise series. 

6. BLOCK CORRELATION 

Two time sequences can be correlated by any of three 
different methods. These will be discussed below. 

The first method blocks the length of each input sequence. 
Correlation of these blocks will provide an output with a 
variance given by Equation 17. The linear dependence of the 
correlation variance on (N- | C | ) is shown in Figure 5 for zero 
mean Gaussian noise. Figure 5 shows the correlated output at 
C =0 (zero time lag) is scaled by (N-|0|)=N. Correlation output 
at a time lag of ±N are scaled by (N-|n|)= 0. The linear 
dependance of the correlation output variance on the time lag 
does not allow a constant detection threshold, which is needed 
to automatically determine the location of the correlation 
peak. 

A second method of correlation blocks a fixed length 
sequence x(n) and correlates it with an arbitrarily long 
sequence y (n) . The correlation output for this algorithm has 




( 17 ) 



where N = number of data points 
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a variance of Na^. Clearly, no scaling of the correlation 
output variance as a function of time lag would occur. 

A third method is scaling applied to the output of the 
first technique. The correlator output is multiplied by a 
weighting function, l/(N-|f|)^^^ to correct for the lag 
dependency of the variance. The selection of a constant 
detection threshold will be addressed in Chapter V. 
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III. FREQUENCY DOMAIN BASED CORRELATION 



A. FOURIER TRANSFORM 

Frequency domain techniques can be used to perform fast 
correlation rather than the time domain techniques previously 
discussed. The Fourier transform translates the discrete time 
sequence x (n) into the frequency domain . 

For a finite duration time series x(n), the discrete 
Fourier transform (DFT) is defined by the pair of equations 
[Ref 5:p. 100] 



X{k) = ^x(n) e ^ , 0 ^ ic <: N-1 

n~0 

( 18 ) 

'I ^ ^ ^ ) 12k 

x{n) = — Y X{k) e ,0 ^ n ^ N-1 . 

^ ^0 



Equation 18 is used to transform the N discrete sample values 
into N frequency domain coefficients. 

The following discussion is based on [Ref. 5:p. 286]. The 
DFT is a computationally intensive algorithm if implemented 
directly. The direct evaluation of the DFT requires N^ complex 
multiplications. The time required for the evaluation of a DFT 
is proportional to N^. 

Algorithms have been designed that take advantage of the 
symmetry and periodicity of the DFT, reducing the amount of 
computation required to solve the DFT. Any of these algorithms 
are referred to as the FFT. 
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Generally, the FFT requires (N/2) log^N complex 
multiplications. This smaller computational workload becomes 
significant as N, the number of data points in the time 
series, becomes larger. 

B. NORMALIZED FREQUENCY DOMAIN CORRELATION 

In some radar receivers the processing of pulse sequences 
occurs in the frequency domain. Therefore, a frequency 
representation of the time domain pulses may already exist. 
The frequency coefficients do not necessarily need to be 
transformed back into the time domain to perform the 
crosscorrelation . 

Frequency domain data can be crosscorrelated by 
multiplying one set of coefficients by the complex conjugate 
of the second set of coefficients. The product is translated 
back to the time domain through the use of the inverse FFT. 
Frequency domain based normalized crosscorrelation is 
described by 
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W-1 J ( -^ ) iJc 

P^(i) = . 0.1. N-1 



\ 



N-1 



N-1 



J2 k(^) P X E I' 

Jc=0 N k^o 



( 19 ) 



= 0 Otherwise 

where X{k) = FFT of the series x{n) 

Y(k) = FFT of the series y{n) 

* = conjugation 

p^(i) = normalized correlation function . 

To avoid circular correlation, both time domain pulse 
sequences are zero padded to Nj data points 



where = number of data points in x{n) (20) 

Ny = number of data points in y(n) . 

To take advantage of the processing speed of FFT algorithms, 
Nj is made a power of two by increasing its length with 
additional zero padding 



2- ^ Ny - 1 ^ 21 ) 

where m is an integer . 

Prior to any zero padding, the dc component is removed to 
create zero mean pulse sequences. Both sequences, each of 
length N, are then transformed to the frequency domain through 
the FFT and processed using the frequency domain correlation 
technique. A MATLAB implementation of frequency domain 
crosscorrelation is given in Appendix B (FreqCorr? .m) . 
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C. TIME DIFFERENCE OF ARR.IVAL (FOURIER DOMA.IN) 

The normalized crosscorrelation function implemented in 
the frequency domain (Ec[uation 19) is used to estimate the 
TDOA. The crosscorrelation of two similar signals will produce 
a dominant peak in the output of the correlator. This peak is 
used as an estimate of the TDOA between the two signals. For 
two signals that have not been corrupted by noise, the 
location in time of this peak is easily determined. Pulsed 
signals that have been distorted by noise produce a noisy 
correlated output. The output noise can mask the TDOA peak in 
low SNR conditions. Consequently, the probability of selecting 
a noise sample instead of the true TDOA peak increases with 
decreasing SNR. 

Figures 7,8 and 9 show the frequency domain based 
crosscorrelation of two pulse sequences as a function of 
decreasing SNR. A reception is simulated with a pulse burst 
demodulated by one receiver fed into one channel, and a pulse 
burst demodulated by a second receiver fed into the second 
channel of the correlator. Pulse burst one has a transmission 
delay of zero. Pulse burst two has a delay of 50 time units. 
Therefore, the total TDOA is 50 time units. The TDOA peak is 
clearly seen at the high SNR of +7dB and +3dB. At OdB SNR, 
noise peaks exceed in amplitude the true TDOA peak causing an 
incorrect estimate if one was selected. 



22 



normalized FFT correlation : channel 1 channel 2 




Figure 7. Normalized frequency domain based crosscorrelation. 
SNR +7dB. 




SNR +3dB. 

Visually estimating the location of the peak with the 
maximum amplitude in the output of the correlator is a easy 
task. The decision rule states that the point with the maximum 
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normalized FFT correlation : channel 1 vs channel 2 




Figure 9. Normalized frequency domain based crosscorrelation. 
SNR OdB. 

amplitude is selected as the TDOA estimate, provided the 
amplitude exceeds a predetermined threshold. 

When many blocks of data must be processed, visual 
inspection of the correlation output is not possible. Chapter 
V discusses the design of a threshold for automated data 
processing. 
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IV. NOISE AND SIGNALS IN RECEIVER 



A. NOISE AND SIGNALS IN THE RADAR DETECTOR 

The purpose of a radar receiver is to process and detect 
the reflected energy from a target being tracked by the radar. 
An I/Q demodulator is assumed in the radar receiver, located 
after the last intermediate frequency (IF) amplifier and prior 
to the video signal display. Data to be processed for TDOA can 
be obtained from three locations on the I/Q demodulator. These 
locations include the envelope output, the envelope squared 
output and the I/Q channel data. 

White Gaussian distributed noise is assumed to be added to 
radar pulses during their transmission through the atmosphere 
and reception by the receiver. This noise contributes to the 
degradation of the SNR and is caused by fluctuations in the 
target radar cross-section and the loss of signal energy due 
to the propagation distance between the target and the 
receiver . 

At typical radar frequencies the front end amplifier of 
the radar receiver contributes to this noise power. The noise 
power for a thermally noise limited receiver is defined as Ng, 
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^ ^sys B 

where k = Boltzmanns constant 
B - IF bandwidth 

( 22 ) 

T = T +T 

sys ar e 

T^j. = receive antenna temperature 
- effective noise temperature 
of the receiver . 

The total noise is modeled as zero mean Gaussian distributed 
noise with a probability density function (pdf) of 



Py-(x) =■ 



2o2 



( 23 ) 



where is the noise variance due to all noise sources (i.e., 
receiver front end^ transmission noise) . 

For an operational radar, either a pulsed signal plus 
Gaussian noise or Gaussian distributed noise alone is assumed 
at the input of the I/Q demodulator. The probability 
distribution of the demodulated signal is dependant on which 
of the three outputs of the I/Q demodulator is used. 

1 . Coherent detection 

If the exact frequency and phase of the received pulse 
burst is known, and matches the frequency and phase of the 
local oscillator in the receiver, then the signal can be 
coherently detected. Coherent detection simplifies the 
probability description of the signal and noise at the output 
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of the I and Q channels. The I and Q sequences have an 
independent Gaussian distribution, and can be described as 

1,0 " ,WQ(n)} ,a^) (24) 

where iUj{n) , niQin) = I and Q mean values respectively 

= noise variance due to all noise sources . 

Both the I or Q channel can be processed using the correlation 
algorithm followed by threshold detection. 

Due to the passive nature of signal reception and 
demodulation, the exact frequency and phase of the received 
pulse burst are not available at the receiver. Therefore, 
coherent detection is not feasible for TDOA estimation. 

2 . Envelope Squared Detection 

For a zero mean Gaussian distributed noise input to 
the I/Q demodulator, and each have a chi-squared 

distribution with one degree of freedom. Their pdf' s can be 
described by 

~ — ^ for y i 0 

tJy/2/ny 

= 0 for y < 0 (25) 

where = variance of I and Q 
y = or 0^ 

The envelope squared output is the sum of and and is 
therefore chi-squared with two degrees of freedom. This pdf 
can be described by 
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foz z i. 0 






2o^ 



where z = + 0^ 



( 26 ) 



Envelope squared detection is attractive because the Hilbert 
transform can be used to enhance the TDOA peak in the 
correlator's output. Unfortunately in a low SNR environment 
the noise power will also increase as the noise contribution 
is squared. 

3 . Envelope Detection 

A common form of receiver demodulation at IF 
frequencies is the envelope detector. For all simulations in 
this thesis, envelope detection will be assumed. The output of 
the envelope detector is the modulation envelope of the 
received pulse burst. For Gaussian distributed noise at the 
input to the I/Q demodulator, the noise at the envelope output 
is Rayleigh distributed with pdf given by 



Py(u) = -^e 20% u > 0 

o2 

= 0, otherwise 

where = variance of Gaussian noise 
u = ^/z = y/j2 + ^ 



( 27 ) 



For a pulsed sinusoid plus Gaussian noise at the input to the 
I/Q demodulator, the pdf of the envelope at the output is 
Rician distributed. The Rician pdf is given by 
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Pair) 




20'^ 



I,{^) , r i 0 
otherwise 



where A = amplitude of received sinusoid 
= variance of Gaussian noise 



( 28 ) 



where ^ gzcos(e) ^ 

® 2n Jo 

and loiz) = modified Bessel function of the 
first kind of zero order . 



(29) 



Figure 10 shows a block diagram of the I/Q demodulator and the 
various pdfs in the radar receiver. 



B. ENVELOPE CORRELATION OF RAYLEIGH NOISE 

Rayleigh distributed noise is produced at the output of 
the envelope detector under noise only conditions. This noise 
is then processed by the correlator during the TDOA estimate. 
In all follow on discussions and simulations, the envelope 
processing scheme is used. 

The autocorrelation of Rayleigh noise as discussed earlier 
produces a large triangular bias in the output which decreases 
as the time lag increases. This bias is produced by the dc 
component of the Rayleigh distributed noise. 

This triangular function contains no information relative 
to the TDOA estimate. The large numerical values in the bias 
can potentially limit the dynamic range in the FFT based 
correlation. A constant, detection threshold cannot be 
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Q Channel 




Figure 10. I/Q demodulator in the radar receiver. 

implemented at the output of the correlator with this 
triangular function present. For these three reasons, it is 
advjmtageous to force the Rayleigh or Rician distribution to 
have a zero mean. This will remove the triangular bias formed 
by the correlation process. This modification of the received 
pulse burst does not degrade the TDOA estimate. For all 
realizations (i.e., blocks of data), the dc component of the 
Rayleigh noise will be removed. 
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1. Shifted Raylttigh Noistt 

Forming a Rayleigh distributed time series x(n) and 
svibtracting its expected value produces a shifted Rayleigh 
sequence. This algorithm is described by 



x^{n) = x( 



n) -o|| 



(30) 



where x^^in) = shifted Rayleigh noise sequence 

o - Gaussian noise standard deviation 



The shifted Rayleigh time series has a zero mean. The new time 
series and the corresponding time domain based autocorrelation 
function are shown in Figure 11. 




(a) (t>) 



Figure 11. (a) Shifted Rayleigh time series, (b) Normalized 

shifted Rayleigh time domain based autocorrelation. 



Histogrzuns of a Rayleigh pdf and a shifted Rayleigh 
pdf are shown in Figure 12. Clearly, the effect of removing 
the dc term from the Rayleigh distributed time series creates 
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Figure 12. (a) Rayleigh pdf. (b) Shifted Rayleigh pdf. 

an autocorrelated output that does not have the limiting 
triangular bias. The impulse at zero time lag dominates the 
correlation output. A constant detection threshold, at least 
over a small number of range bins, could be implemented with 
this type of correlation output. 

2 . Frequency Domain Based Correlation of Shifted Rayleigh 
Noise 

The frequency domain based autocorrelation of Xi(n) is 
now examined. Sequence x^ (n) is zero padded as described in 
Equation 21. A frequency domain series is created using the 
FFT as described by Equation 18. The autocorrelation estimate 
is derived using Equation 19 and is plotted in Figure 13. The 
results of frequency domain based autocorrelation are very 
similar to the results of time domain based autocorrelation 
using shifted Rayleigh noise. 
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normalized shifted Rayleigh spectral autocorrelation 




Figure 13. Normalized frequency domain based autocorrelation, 
3 . Spectral Interpolation 

With the incorporation of the FFT in modern radar 
receivers, it is possible that the received pulse burst has 
been transformed into frequency data as a byproduct of 
detection. Processing time would be lost in converting this 
data back into a time series to be correlated. Of course, 
there is no guarantee that the dc component has been removed. 
In this section an algorithm is given to perform frequency 
domain based correlation on a Rayleigh distributed noise 
sequence while avoiding the triangular correlation bias 
created by the (possible) dc component. 

A Rayleigh distributed noise series x (n) , N points 
long, is transformed into an ordered sequence of N complex 
coefficients using the FFT algorithm. In general, a sequence 
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X(k) of coefficients is indexed from zero to N-1 and is 
described by 

X{k) = FFT { xin) } , 0 ^ k iiV-1 

where k = frequency index . 

The FFT creates the zero indexed term X(0) by summing all the 
elements in x(n) . X(0) is real and can be viewed as the dc 
component of x (n) . By removing the X(0) coefficient from the 
frequency sequence we have in effect removed the dc bias from 
the time series x (n) . At this point, the modified Fourier 
transform X (k) approximates the DFT of a shifted Rayleigh 
series . 

Time domain correlation of an N point series produces 
2N data points. Frequency domain based correlation of an N 
point sequence produces N data points. To force frequency 
domain based correlation to equal time domain based 
correlation (i.e., avoid circular correlation), the number of 
terms in the frequency series must be doubled. The frequency 
series is expanded by storing each pair of complex frequency 
coefficients at a location twice the value of the original 
index . 

As a first approximation, the data points between 
adjacent coefficients are linearly interpolated in the 
expanded array. Figure 14 illustrates the effect of FFT 
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interpolation and the zeroing of the dc terms. Correlation 
results as described by Equation 19 are shown in Figure 15. 



Re{ X(k) } or lm{ X(K) } 



A 



X(-N/2) X(-2) 1 

1 1 1 


1 x(- 1 ) 

ill. 


X(1) 1 X(2) X(N/2) 

i 1 1 1. 1 ^ 


H 1 1 


I 1 — t — ^ 


T 1 


1 1 1 



-N/2 -4 -3 -2 -1 0 1 2 3 4 N/2 



I = interpolated coefficient 
X(k) = Fourier coefficient 



Figure 14. FFT interpolation, and zeroing of the dc 
coefficient . 




Figure 15. Normalized frequency autocorrelation through 
interpolation . 
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Clearly, the triangular bias in the correlation output has 
been removed. A small impulse at a time lag of -N is observed, 
which can be disregarded. 

For block processing, typically only correlation 
products within ±10% of the zero time lag position are 
considered to be statistically reliable. Therefore, the TDOA 
estimate is formed within this correlation band. The impulse 
at -N falls outside of the TDOA evaluation area and will not 
cause a false detection. 

C. CORRELATION OF PULSES IN ADDITIVE RAYLEIGH NOISE 

Time domain based correlation of a pulse train in additive 
Rayleigh noise should produce the same result as frequency 
domain based interpolation followed by an equivalent 
correlation. A comparison was performed using a 50% duty cycle 
pulse train in additive shifted Rayleigh noise. 

Time correlation was performed on the pulse burst. The dc 
component of the noisy pulse train is removed forming a zero 
mean signal. The time domain based autocorrelation of this 
sequence is shown in Figure 16 (a) . 

The pulse burst was also transformed into the frequency 
domain using the FFT with the dc component intact. N Fourier 
coefficients are expanded and the X(0), X(±l) coefficients 
zeroed out. The sequence is correlated after interpolation 
using Equation 19 and plotted in Figure 16(b). 
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Figure 16. (a) Normalized time domain based autocorrelation, 
(b) Normalized frequency domain based interpolation and 
autocorrelation. Zero time lag between each channel. 



Clearly, the two correlation algorithms produce similar 
outputs. Both have a dominant peak in their output for TDOA 
estimation. Frequency domain based correlation displays a 
exponential decay of the sidelobes about the main peak. This 
effect is not observed in time autocorrelation, and is caused 
by the first order interpolation used in the frequency domain 
based method. A more robust interpolation algorithm will 
eliminate the exponential decay, yielding the desired linear 
decay . 

Next, time and frequency domain based crosscorrelation are 
performed on two pulse sequences. One pulse sequence lags the 
other by 20 units. Both time and frequency domain based 
crosscorrelation of the two sequences are plotted in Figure 
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17. Clearly, both correlation algorithms produce a peak at 20 
units. The nonlinear decay of the sidelobes surrounding the 
dominant peak is again seen in frequency domain correlation. 





(a) (b) 



Figure 17. (a) Normalized time domain based crosscorrelation, 
(b) Normalized frequency domain based interpolation and 
crosscorrelation. 20 unit time lag between each channel. 



In summary, the frequency domain based interpolation and 
crosscorrelation algorithm produces an output comparable to 
time domain based crosscorrelation. It offers a method to 
correlate pulse bursts collected in the frequency domain. 
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V. THRESHOLD DESIGN FOR ENVELOPE CORRELATION 



A. INTRODUCTION 

A threshold at the output of a correlator is used to 
estimate the location of a TDOA peak in a constant false alarm 
rate (CFAR) receiver. The value of the threshold is a function 
of the correlation output mean and variance. 

B. CORRELATION OUTPUT AND CONSTANT VARIANCE 

If the variance of the correlator, when driven by two 
signals at the input, maintains a fixed value about some mean, 
a constant threshold above the mean can be used. This allows 
the design of a constant false alarm rate detector. The 
largest correlation peak that also crosses the threshold is 
defined as the TDOA estimate. The constant threshold design is 
the easiest to implement. 

C. CORRELATION OUTPUT AND TRIANGULARLY SHAPED VARIANCE 

The crosscorrelation of two zero mean noise sequences with 
the same number of data points, produces an output whose 
variance has triangular form. Hence, a constant detection 
threshold is not easily implemented. Instead, a threshold must 
be designed that has a slope that matches the slope of the 
correlation output. The design of a sloping threshold is 
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difficult, particularly if the slope of the correlation peak 
changes with time (i.e., changing noise environment). 

A solution is to mathematically force the triangular 
correlation variance to be flat. This is done by using a 
bowtie correction of the form 1/ (N- | C | A bowtie correction 
normalizes the correlation output to have a constant (lag 
independent) variance. This design suffers from time changes 
in the correlation output forcing the bowtie correction to be 
recalculated. Note that the amplitude of the TDOA peak will 
itself be dependant on the type of correction to the 
correlation output. 

D. THRESHOLD BASED ON ZERO LAG 

The zero lag threshold method uses the variance of the 
zero lag correlation output to calculate a constant detection 
threshold. This is indicated in Figure 18. 

The correlation output becomes statistically less reliable 
the further removed the correlation products are from the zero 
time lag position. To use correlation intelligently, the 
correlation output is typically examined ±10% of the distance 
from the zero time lag position. Consequently, only 
correlation output ±0 . IN from zero time lag will be used in 
TDOA estimate simulations. 

The detection method in this thesis uses a constant 
threshold against the correlation output. The fixed threshold 
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was selected because the correlation output variance is, to a 
first approximation, nominally flat within the ±0.1N boundary. 

The detection scenario that drives this thesis defines a 
minimum radar pulse width of one microsecond (^.sec) and a 
burst length of four millisecond (msec) . A 50% duty cycle 
radar pulse is assumed. There are 2000 pulses per burst, which 
when correlated, will generate an output that is SOOO^lsec 
(±4000|lsec) in length. If ±10% of the correlation output 
around zero time lag are compared to the threshold value, 
±400jlsec will be tested. For radar pulses traveling at the 
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speed of light, a maximum distance of ±120km is therefore 
examined to form a TDOA estimate. 

Correlation values at lags -O.IN to 0 . IN will be compared 
in magnitude to the zero time lag based calculated threshold 
value . 



E. THRESHOLD CALCULATION FOR CFAR 

The detection threshold is calculated using the output 
terms of the crosscorrelation algorithm when noise only is 
injected into the correlator. For convenience, Equation 8 is 
repeated here, 

AT-UI-l 

= E yi^i 

(32) 

and are zero mean Rayleigh {shifted) 

distributed noise sequences . 

The correlation function at all lags of interest is thought to 
be Gaussian distributed. According to Equation 32, every 
estimate is the sum of a fairly large number of products 
(i.e., central limit theorem) . The average value of the 
correlation output between the ±0.1N endpoints is 



t 



1 

2a + 1 






E 

7 = -x 



r^yU) 



where a = number of data points between the 
0 and O.IN lag positions In the 
correlation output . 



( 33 ) 



The variance of r,y(l) is calculated from 
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2a + 1 



(34) 
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(r,^(J) 



t )2 



For Gaussian distributed noise the probability of false alarm 
(Pfa) is a function of the detection threshold 



p = 1 r 



20fc 



dx 



(35) 



where T = detection threshold 



For a given Pf^, the threshold can be calculated from 

T = erfc-^ (Pfa) Ot 



(36) 

where erfc = complementary error function . 

Equation 36 relates the threshold to the Gaussian distributed 
noise variance. 



F. VERIFICATION OF CFAR THRESHOLD AND ROC CURVE 

Using the above ecpaations, thresholds were calculated for 
the time correlation detector given in Equation 32 for three 
different Pf,'s. A MATLAB simulation was performed to compare 
the measured false alarm rate of the time correlation detector 
to the theoretical false alarm rate. 

1 . Pf. Simulation 

Two 100 point sequences of zero mean Rayleigh noise 
were injected into the correlation detector. The output of the 
correlation detector between the -O.IN and +0 . IN lag points 
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was compared to the calculated threshold to measure the false 
alarm rate. The number of points in the correlation detector 
output that exceeded the threshold for a given Pf. were 
recorded. The experimental Pf^ is given as the ratio of 
improper threshold crossings to the total number of monitored 
range cells of the correlator output (see Table 1) . For each 
Pfa, 500 realizations were performed. Both theoretical and 
experimental results are listed. 



TABLE 1 

THEORETICAL PROBABILITY OF FALSE ALARM 


Simulation # 


» 

11 

O 

o 


Pj.=0.001 


Pf,=0. 0001 


1 


0.0100 


0.00076 


0 


2 


0 . 0099 


0.00090 


0 


3 


0.0110 


0.00085 


0 



The data in Table 1 show that a detection threshold 
can be established for a correlation detector, using the 
assumption of a Gaussian distributed output for Pf,' s of 0.01 
and 0.001. For the lower P^^ value of 0.0001, the threshold 
calculation fails, but in a positive sense. The measured Pf^ 
of 0 is below the theoretical Pf^ of 0.0001. Either the 
correlation variance is smaller than Equation 36 predicts, or 
more terms must be summed in the correlation output to better 
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approximate the Gaussian pdf as indicated by the central limit 
theorem. 

2 . TDOA Simulation 

The output of the envelope detector will either be 
Rayleigh or Rician (non-central Rayleigh) distributed. In 
those portions of the data where there is no signal energy, 
the output will be Rayleigh distributed. Where the signal is 
present, the envelope detector output will be Rician 
distributed. 

The purpose of the time domain based correlation 
detector is to measure the TDOA of pulsed sequences embedded 
in Rayleigh noise. Noise in low SNR environments will mask the 
TDOA peak. A MATLAB simulation was performed to examine the 
performance of the correlation detector in a decreasing SNR 
environment . 

Two pulse bursts embedded in additive Rayleigh noise 
were created. The second pulse burst lagged the first by three 
time units. Each pulse burst is a sequence of ten pulses, each 
having a period of ten units. Each pulse in the burst had a 
50% duty cycle. 

The SNR established the amplitude relationship between 
the noise and signal sequences. First, two sequences of zero 
mean unit variance Gaussian distributed noise were created. 
Rayleigh noise was generated from the two Gaussian sequences. 
The Rayleigh noise had a mean of 1.2533 and a variance of 
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0.4292. The theoretical second moment of the Rayleigh noise is 
two and was verified in the MATLAB simulation. The second 
moment is defined as the average power of the Rayleigh noise. 
The amplitude of the nonzero elements in each pulse burst was 
then scaled to meet the desired SNR, using the measured 
average power of the Rayleigh noise. 

A detection threshold was calculated based on shifted 
Rayleigh noise using Equation 36 where <5^ was assumed to be 
known (i.e., simulation information). In practical 
applications, <7^ would be estimated from the data. The SNR of 
the two bursts were varied from OdB to 18dB in the different 
simulations. For each SNR, the two sequences were correlated 
using Equation 32. 

Graphs of the correlation output for a SNR of IdB are 
given in Appendix C for Pfa' 3 of 0.01, 0.001, 0.0001 and 
0.00001. The threshold, which is a function of the Pfa and 
noise variance, is shown to increase with decreasing Pfa« 

The correlation peak at a time lag of three was the 
correct TDOA estimate. If the maximum correlation peak in the 
output was the third lag point, a correct TDOA estimate was 
obtained. For each SNR, 30 TDOA estimates were made. The 
correlation function output 100 data (lag) points per TDOA 
estimate. The TDOA estimate was made on ±10 data points on 
each side of the zero time lag position in the correlation 
output. The percent correct TDOA estimate for each SNR was 
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calculated and is presented as a receiver operating 
characteristic (ROC) curve in Figure 19. 

The ROC curve shows that as the SNR increases, the 
performance of the correlation detector increases. The 
performance of the correlation detector in Figure 19 can be 
improved if multiple sequential looks at the same SNR were 
allowed. 
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Figur* 19. Percentage of estimated TDOA correct versus 
SNR for time domain based correlation detection. 
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VI. FFT ERROR ANALYSIS 



A . INTRODUCT I ON 

This chapter discusses several sources of processing 
errors occurring within the FFT, and the effect of these 
errors on the FFT output signal-to-noise ratio. 

The computation of a FFT using fixed point arithmetic, 
generates internal processing errors. These errors are due to 
the fixed number of bits in the trigonometric look up table, 
the truncation of arithmetic results during addition and the 
scaling of results during the individual butterfly 
calculations . 

The error sources can be modeled as independent additive 
noise components at the output of the FFT. Their net effect is 
to reduce the ideal SNR at the output of the FFT. 

Each of these processing error sources will be described 
below. The material in this chapter borrows from [Ref. 7] . 

B. COSINS/SINS TABLE NOISE 

The computation of a butterfly algorithm in the FFT 
involves multiplication by the complex coefficient (twiddle 
factor) 
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(37) 



e 




cos(-^n?) - j sin(-^/n) 

N N 

W” . 



A FFT butterfly is shown in Figure 20. 




These trigonometric coefficients are precalculated, quantized 
to (B) bits (where B is the size of the computer word 
including sign bit) , and stored in memory for future table 
lookup during the FFT processing. The coefficients W® * 1 and 
= ~j cam be obtained in an exact form and therefore do not 
contribute to the quantization noise. Quantization of the 
other (twiddle) coefficients stored in the sine and cosine 
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table introduce noise into the output of the FFT. The noise 
power created by truncation of the sine and cosine 
coefficients is defined by [Ref. 7] to be 




(38) 

where ol = trigonometric noise power 

B = number bits in computer word . 



C. TRUNCATION NOISE 

The DFT of a finite duration sequence (x(n) } is defined in 
Equation 18 and is repeated here for convenience 

-illn 

X{k) = ^ x(n) e ^ , 0 ^ k ^ N-1 . (39) 

The product formed requires four real multiplications because 
both the exponential and the input data are complex numbers, 
B bits in length. Each multiplication can produce a result of 
2B bits which must be truncated from 2B to B bits, hence there 
are four quantization errors for each complex valued 
multiplication. The four quantization errors can be described 
as four independent noise sources. The truncation noise power 
is defined by [Ref. 7] as 
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( 40 ) 






where o\ - truncation noise power 

B = number of bits in product after truncation 



D. SCALING NOISE 

The FFT as an algorithm processes a data vector of length 
N by successive passes at the vector. During each pass, the 
algorithm performs N/2 butterfly operations. Each butterfly 
retrieves two complex numbers from memory, performs the 
butterfly computation and returns two complex numbers to the 
same memory address . The complex numbers returned to memory by 
the butterfly can have a greater magnitude than the two 
complex numbers initially fetched from memory by the 
butterfly. In order for the results of the butterfly operation 
to fit in the fixed word length of the memory of the computer, 
the results must be truncated (scaled) . Scaling is performed 
by dividing by two the entire data array, and is implemented 
in software by a right shift and an increment of the arrays 
exponent register. 

During the (m+1) pass of the data, the butterfly algorithm 
selects two data points A(m) and B (m) and returns to memory 
A(m+1) and B(m+1). The truncation error results from the 
addition of two numbers of like sign in the upper part of the 
butterfly algorithm or the subtraction of two numbers of 
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opposite sign in the lower part of the butterfly algorithm. 
The input to the butterfly is fixed to have maximum real and 
imaginary values of ±1. The output of the butterfly is fixed 
to have a maximum value of ±2, double the input value. To 
prevent overflow of the butterfly output, either a prescale by 
1/2 is required prior to entering the butterfly or an extra 
bit must be available in memory to accommodate the possible 
gain of two. 

One method, automatic prescaling, truncates the least 
significant bit of the butterfly output and represents 
processing noise. In the case for which no scaling is 
necessary for the pass, this truncation represents significant 
processing error (noise) . 

A second method, data dependent scaling, is performed only 
if a butterfly in the data pass overflows without scaling. If 
any word in memory has an overflow bit set, all words are 
shifted right as they are fetched from memory for the next 
pass of the FFT. The total number of right shifts, p, executed 
during the FFT process is used at the output of the FFT as a 
scale factor of 2 ^. This scale factor is used at the output of 
the FFT so that the total processing gain is maintained. 
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The noise power created by the scaling of the data is 



defined by [Ref. 7] 



,2 _ 






where o| = scaling noise power 



( 41 ) 



E. OUTPUT SIGNAL TO NOISE RATIO 

Figure 21 shows an error model of a butterfly in an FFT. 



q^(m) 




Figure 21 . FFT error model . 
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The input data terms to a butterfly can be described by 



A (m) = A im) + Nj^ (m) 

B{m) = B{w) + Ng (m) 

( 42 ) 

where Nj^, Ng = additive noise terms 
A{m) , B{m) - input data terms 
A(,m) , B{m) = noise corrupted data . 

without the effects of noise the butterfly algorithm produces 

two outputs described by 

A{m+1) = A(m) + W(m)B(m) 

B{m+1) - A{m) - W{m)B(m) ( 43 ) 

where W{m) - twiddle factor . 



The error due to the A component will be calculated. The error 
due to the B component is statistically equivalent to the 
error in the A component. Including the effects of noise, the 
output of the butterfly can be described by 



A(/77+1) = [A(m) + Nj^im) + (7g(/n)] + [B{m) + Ng{m) ^-q^im)]- 

[W(m) + g„{m)] + Qj-im) 



where dgim) = scale noise 

g„{m) = sine/ cosine table noise 
gj.{m) = truncation noise (in multiply) 



( 44 ) 



The total power at the output of the butterfly is then 



55 



( 45 ) 



£■{^(^+1) A’{m+1)} = E{A{m)A‘{m)} * E{B{m) B' (m) } + 

E{Nj^{m) (m) } ^ E{Ng(m) Nb (m) ] + 
E{gs(m) Qs (m) } E{qg(m) (m) ] + 
E{B{m) B' (m) } Eig^^im) g„[m) } + 
E{gj.(m) gj (m) } . 



Redefining the terms of the total power leads to 



E[A(m+l)A’ (m+1) } 

where S'^ {m) 
ol (m) 
a^{m) 
al(m) 



S^(m) + (m) + Nl(m) + {m) + 

20s(m) + (m) al/im) + Ojim) 

input signal power 
scaling noise power 
truncation noise power 
sine/ cosine table noise power 



Combining like terms results in 

£’{A(m+l)A*(m+l) } =2S^(m) ^ 2Nl{m) +2a^(/n) + 

(m) oj,(m) + Ox(m) . 

This equation describes the total power as a function of the 
signal power and the individual independent uncorrelated noise 
terms . The total power output from the butterfly can also be 
described as a function of its input signal power and the 
total additive noise power generated within the butterfly as 

E{A{m+l) A’ {m+D ] = E{A{m+i) A’ (m+l) ] + 

E{Nj^{m-^l) Na (/n+l) } 

Redefining these terms 

E{A{m+l)A'(m+l)} = S^{m+1) + Nlim+l) . 



( 48 ) 



( 49 ) 
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Equation 47 and 49 are equivalent. The signal and noise power 
for the (m+1) pass can be formed from those terms in the 
proceeding pass 



+ iv|(/77+l) = 2S^ irn) + (m) + 2o|(/n) + 

im) olim) + Or(/n) . 



(50) 



Equating the signal terms in Equation 50 yields 

(m+1) = 2S^ (m) . (51) 



Equating the noise terms in Equation 50 yields 

N^im+l) = 2N^(m) + 2o|(/7?) + (m) a^(m) + Q%(m) . (52) 

Equation 51 can be rewritten to show how the input signal 
power increases as a function of m as it passes through the m'^*' 
stage of a butterfly of an FFT 



S2 (m) = 2'"S2, 

here we modeled 5^(0) =5^. 



(53) 



The signal power input to the butterfly is defined as 

S2 = l2-2(s - *> . (54) 

9 

Rewriting Equation 52 to include this exponential increase of 
S* (m) gives 
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wj(/n+l) = 2Nl{m) + 2al(m) + 2‘”S^alf{m) + a^im) . (55) 

Equation 55 shows the relationship between the noise power of 
a butterfly stage as a function of the noise power from the 
proceeding butterfly stage. The noise and signal power can now 
be calculated using Equations 53 and 55 for m passes through 
a butterfly. Once calculated^ the signal to noise ratio can be 
formed. 

F. RIGHT JUSTIFIED DATA 

The length of a word in the computer is defined to be 16 
bits (B) and accepts a 12 bit (b) data word. Right justified 
data places the least significant bit (LSB) of the b bit data 
word in the LSB of the B bit processor word. The output noise 
to signal ratio is now calculated for ten passes through a 
butterfly using right justified data. The following 
assumptions are used in the calculation; 

1. The input signal is incoherent. 

2. Scaling is performed every other pass after the most 
significant bit (MSB) of data reaches the next to MSB 
of the processor. 

3. The first two passes do not use the trigonometric 

tcdDle . 

After the tenth pass (i.e.^ FFT size of 1024), it is shown 
(Appendix D) that the noise to signal power is 
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(56) 





? 1 O x 

40 ^ + + 

4 52 



1 



a 



2 

s 



51 52 ’ 



B and b are previously assumed to be 16 and 12 bits in length 
respectively. Substituting these values into Equation 56 




3 ^ 

_l2"8 3 4 ^2'® 1 2~® 

9 9 9 



. (57) 



The noise to signal ratio reduces to 




-67.5 db - 95.1 dS -71.1 dS -79.1 dB (58) 

input trig truncation scaling 

noise noise noise noise . 



Clearly, the dominant degradation noise is the truncation 
table noise (-71.1dB) , followed by the scaling noise to signal 
ratio (-79.1dB) . For right justified data, the signal to noise 
ratio is calculated to be 65.71dB. 



G. LEFT JUSTIFIED DATA 

Left justified data places the most significant bit (MSB) 
of the b bit data word in the MSB of the B bit FFT processor 
word. The output noise to signal ratio is now calculated for 
10 passes through a butterfly for left justified data. The 
following assumptions are made 
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1. The signal is incoherent. 

2. Data is left justified on input and scaling is 
performed every other pass starting with the first 
pass . 

3. The first two data passes do not use the 
trigonometric table. 

It can be shown that after the tenth pass through the 
butterfly that the noise to signal ratio is 





+ 4o^ + 



1 



O 



4 52 




( 59 ) 



B and b were previously assumed to be 16 and 12 bits 
respectively. Substituting these values into Equation 59 we 
obtain 





1 J^2-32 
4 12 





( 60 ) 



The noise to signal ratio reduces to 




-67.5 dS -95.1 dS -89 . 1 dB -77 . 1 dS ( 61 ) 

input trig truncation scaling 

noise noise noise noise . 



Clearly, the worst noise contribution to the SNR is the 
scaling noise (-77.1dB), followed by the truncation noise (- 
89.1dB) . For left justified data, the signal to noise ratio is 
calculated to be 67.0dB. 
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Further analysis shows that for large FFT's (i.e., in the 
range of 2^® and above) , the trigonometric table noise is the 
dominant noise source. A simulation was performed using an 
AMDAHL machine to exaunine the error generated in performing a 
FFT and an iFFT. Figure 22 shows the experimental setup. 



AMDAHL Machine 




Figur* 22. FFT processing error flowchart [Ref. 8]. 

The error that was generated was taJDulated. Figure 23 
illustrates the error as a function of FFT length. 

H. SUMMARY 

For either left or right justified data, the dominant 
noise at the output of the 1024 point FFT is the scaling or 
truncation noise. The SNR at the output of the FFT, for left 
justified data is 67.0dB and exceeds the SNR for right 
justified data which is 65.71dB. Clearly, left justified data 
maximizes the computational SNR at the output of the FFT for 
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MAXIMUM ERROR IN SINGLE 
PRECISION FFT 




FFT LENGTH 



Figura 23. FFT error as a function of transform length 
[Ref. 8] . 
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the values of b, B and N used in the above examples. 

For a right justified data set transform size of N=2^^ 
(8192 points) , the noise to signal ratio after the 13'"*' pass 
can be written 



+ (2^5 +2“ + 2 ^ 2 ) _ 52 <j 2 ^ 2047Or + 170a 



(13) 



21^52 



(62) 



This equation reduces to 






(13) = 



7 107 * 

+ 5 . 50 ^ + d:_I 
52 4 52 



48.2 52 



(63) 



Comparing Equation 56 to Equation 63, we see that by 
increasing the transform size from 2 ^° to 2^^ increases the 
noise to signal ratio. The trigonometric table noise is 
increasing faster than either the scaling or truncation noise. 
The signal to noise ratio is reduced, from 65.71dB to 65.68dB. 
As the transform size is increased to accommodate larger data 
sets, the trigonometric table noise becomes the dominant error 
term. 
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VII. SUMMARY 



A. CONCLUSIONS 

A frequency domain based algorithm was developed and 
tested to estimate the differential arrival time of a pulsed 
radar signal collected by two passive sensors. For 
convenience, the actual implementation was performed through 
time domain processing even though frequency domain processing 
is advocated. The performance of the algorithm is 
characterized by a ROC curve as a function of SNR. For 3dB SNR 
and a of 0.01, the probability of making a correct TDOA 
estimate exceeds that of an incorrect estimate for a single 
look. At 18dB SNR, the probability of making a correct TDOA 
estimate is 100 percent for all Pf.'S. If multiple looks (i.e., 
multiple bursts) are allowed, the probability of making a 
correct TDOA decision at each SNR will increase. 

An I/Q demodulator is assumed in the radar receiver. The 
pdf of the signal driving the correlation detector is 
determined from where it is obtained in the I/Q receiver. The 
pdf, depending on selection, will be zero mean Gaussian versus 
non-zero mean Gaussian, or chi-squared versus non-central chi- 
squared, or Rayleigh versus Rician (non-central Rayleigh) . 
This thesis assumes an envelope detector at the output of the 
I/Q demodulator. The envelope detector output has a Rayleigh 
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or Rician pdf and drives the correlation detector. The 
calculation of a CFAR threshold assumes a Gaussian pdf at the 
output of the correlation detector. For small time lags, 
summation of the output of the correlator may produce a 
probability distribution that deviates from Gaussian. This 
deviation would produce a biased threshold calculation. To 
minimize the bias, a sufficient number of terms must be summed 
at the output of the correlator to allow the Gaussian 
approximation . 

If the received pulses are collected in the frequency 
domain, a spectral domain based correlation algorithm can be 
implemented. An algorithm is given in this thesis. 

For any digital signal processing algorithm that uses the 
FFT, processing errors must be considered. For large transform 
sizes the trigonometric noise power dominates. The length of 
the trigonometric coefficient word affects the degradation of 
the SNR at the output of the FFT. The larger this word the 
smaller the noise power. For a given transform size, left 
justified data will have a higher output SNR than right 
justified data. 

For a correct TDOA estimate, the correlation algorithm 
requires the reception of the leading pulses in the radar 
pulse burst. If the pulses are received in an adequate SNR, 
the correlation algorithm will produce a well defined peak 
allowing the TDOA estimate. Should the environment degrade the 
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destructive interference due to 



received signal (i.e., 
multipathing^ weather or terrain ) , the correlation algorithm 
should be discarded in favor of the traditional angle of 
arrival (AOA) algorithm (i.e.^ maximize received energy). 

If the radar uses a staggered PRF^ the reception of the 
first pulse or several pulses has little impact on the 
position of the correlation peak. Under this condition^ the 
TDOA algorithm is superior to the AOA algorithm, 

B . RECOMMENDATIONS 

A ROC curve should be constructed for the correlation 
algorithm used in this thesis for multiple TDOA looks. An 
intensity display could be designed that would plot as a 
function of time (i.e,^ snapshots)^ those lag points chosen as 
TDOA estimates. Patterns displayed on the intensity plots 
would allow the user to visually determine the correct TDOA 
(i.e.^ incoherent averaging). 

The performance of the spectral domain based correlation 
algorithm developed in this thesis must be further quantified. 
Sets of ROC curves should be obtained for different SNR' s 
(i.e,^ SNR variations between channels). 

The simulated FFT error graph should be validated by 
measuring the error performance of a dedicated properly 
dimensioned FFT, 
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APPENDIX A. CORRELATION MEAN AND VARIANCE 



Two zero mean, independent time series are the inputs to 
a correlator. This appendix will derive the expected value and 
variance of the crosscorrelation function output. In this 
thesis, the noise is zero mean (i.e., shifted) Rayleigh noise. 

The crosscorrelation function r^y(l) is defined in Equation 
8 and is repeated here for convenience 



A. EXPECTED VALUE 

Xi and Yi are assumed to be independent, zero mean 
sequences 




x(i) y{i + 1) 



( 64 ) 



i=0 



E[x^} = E{yi) = 0. 



( 65 ) 



The expected value of r,y(l) is 






W-I-UI 



E{r^u)} = E{ = E 



i=0 




( 66 ) 



= 0 
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B . VARIANCE 



The random variable z has a variance defined to be 

ol = £{ (z - E{z} ) 2} 

because E{z} =0 for a zero wean sequence, (67) 

o| = E[z^} . 

The variance of the crosscorrelation function r^y(l) is then 

N-l-|i| N-l-UI N-l-\l\ 

i-0 i=0 i=0 ^ggj 

N-l-\l\ W-l-UI 

= E E • 

i»o j*o 

and are two groups that are independent of each 

other. The expectation operation can be applied to each group 
individually . 

= E E 

i=o J*o 

The two indices i and j define a matrix of values. Two 
contributions have to be considered. 

Contribution 1. When i=j the summation occurs along the 
diagonal of the matrix. Only one summation is used and 
essentially the terms are being squared. These squared terms 
are not independent. To simplify the computation define m 
equal to i and j 
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(70) 



q = ^ E{xl}E{yl.i} 



m-0 



Because and y„ are zero mean sequences, 



E{xl) = ol, E{yl,] = of 



and = Oy = 



Substituting into Equation 70 



W-l-UI N-l-UI 

Cl = ola^ = ^ 1 



w-Q 



171=0 



= o* [ (N-l-|i|) +1] 



= oUN-llj) 



(71) 



Contribution 2. When i^j all the terms except those on the 
diagonal are summed. These terms are independent. For i?tj 



w-i-Ul w-i-Ul 

Cz = ^ E{x^} EiXj} E{y^^^} Elyj^^'^ 

1=0 j=o 

= 0, using Equation 6 5 



(72) 



Therefore, Equation 69 becomes ^ = q + C 



^ ^ (73) 

= <J*{N - \1\) . 
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C. CROSSCORRELATION STATISTICS 



Each data point output by the correlation function 
represents a summation of terms. Using the central limit 
theorem, the output of the correlation function is assumed to 
have a Gaussian distribution. 

In conclusion ^xy(l) - N(0, o‘‘(iV-|i|) 
where N = Normal (Gaussian) distribution . 
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APPENDIX B. MATLAB CODE 



A. USERS GUIDE 

The frequency domain based correlation software given in 
this appendix is written in Pro-Matlab (version 3.5h). All 
simulations were run on a Naval Postgraduate School Sun work 
station using a UNIX operating system. 



B. CORRELATION MATLAB CODE 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Prograun Name : FreqCorr7.m 12 May 91 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear 

clg 

clc 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% PULSE CORRELATION 

% Version 7.0 

% 

% 1. 1 Channel vs Reference Channel 
% or 

% 1 Channel vs Second Channel. 

% 

% 2. Rician Signal and Rayleigh Noise pdf's . 

% (Low SNR signals have power approximately 2.) 

% 

% 3. 0 dc Correlation. Subtract dc from both signals. 

% 

% 4. Normalize correlation output in Fourier domain. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% This MatLab program will either : 

% 1. Correlate a received pulse sequence against a reference 
% pulse sequence. The reference sequence parameters are 
% specified by the user and are used to construct the 
% sequence . 

% or 

% 
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% 2. Correlate two received pulse sequences against each 
% other. 

% 

% 

% Correlation of the pulse sequences is accomplished by 
% multipling in the frequency domain one spectrum against 
% the conjugate of the other spectrum. An inverse FFT is 
% performed on the result to yield the time difference 
% (TDOA) between the two sequences. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



clc 

echo on 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% SELECT CORRELATION OPTION 

% 

% Select either 1 or 2 : 

% 

% 1. Time delay in one receiver channel measured against a 

% reference signal. 

% 

% 2. Differential time delay between two channels. 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
echo off 

option = input (' Select correlation option 1 or 2 ; '); 

clc 



echo on 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Option 1 : Time Delay In One Channel 

% 

% Reference Pulse Train Parcimeter Selection. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
echo off 
if option==l 

wi=input ( 'Pulse Width (seconds) : '); 

pr=input ( 'Pulse Period (T seconds) : ')/ 

nu=input ( 'Number of pulses in pulse train : ')/ 

NumberPointsRef=pr*nu; 

end 
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clc 

echo on 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% 

% Delayed Pulse Train Parameter Selection 

% 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
echo off 

dly=input Channel 1 delay tau <sec.) of pulse train : '); 

wi2=input <' Pulse Width (seconds) : '); 

pr2=input (' Pulse Period (T seconds) : '); 

nu2=input (' Number of pulses in Delayed Pulse Train : ')/ 

NumberPointsDel= (pr2*nu2) +dly / %calc. nmbr pts sig + delay 

clc 

if option==2 

dlyCh2=input (' Channel 2 delay tau (sec.) pulse train : ')/ 

wiCh2=input (' Pulse Width (seconds) : '); 

prCh2=input (' Pulse Period (T seconds) : '); 

nuCh2 = input (' Number of pulses in Delayed Pulse Train : '); 

NumberPointsDel2= (prCh2*nuCh2) +dlyCh2; %nmbr pts sig + delay 
end 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Zero pad both sequences to N3=N2+N1-1. Must then make N3 
% meet the simple formula N3 = 2''m to allow speedy 
% computation of the FFT. 

% WARNING : IBM 80286 MATLAB WILL NOT ALLOW VECTORS GREATER 
% THAN 4000. So an input pulse train of 128 will exceed the 
% system. Solution is to use the SUN workstation or 386/486. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if option==l 

ZeroPadPoints=NumberPointsRef+NumberPointsDel-l / %N3=N2+N1-1 
for m=2:l:19/ %2^ (19) = 524288 point FFT max 

if 2''m >= ZeroPadPoints , 

ZeroPadPoints=2''m; 

break; %when find a 2''m power, exit loop 

end 

end 

end 
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if option==2 

ZeroPadPoints=NumberPointsDel2+NumberPointsDel- 1 ; %N3=N2+Nl- 1 
for m=2:l:19; %2''(19) = 524288 point FFT max 

if 2''m >= ZeroPadPoints , 

ZeroPadPoint 3=2 ''m; 

break; %when find a 2''m power, exit loop 

end 

end 

end 



%%%%%%%%%%%%%%% Create Reference Pulse Train %%%%%%%%%%%%%% 
if option==l 

ReferencePulse=zeros (1 : ZeroPadPoints) ; %zero pad past pr*nu 
for PulseCounter=l : pr ; NumberPointsRef %steps of Ref. period 
for CutUpPulse=0 : pr ; 
if CutUpPulse<=wi 

ReferencePulse (PulseCounter+CutUpPulse) =1.0; 

end 

end 

end 

Refdc=mean (ReferencePulse (1 : NumberPointsRef) ) ; 

ReferencePulse (1 : NumberPointsRef ) =ReferencePulse (1 :NumberPoi 

ntsRef) -Refdc; 

end 



%%%%%%% Create Channel 1 Delayed Pulse Train %%%%%%%%%%%%%% 
DelayedPulseTrain=zeros (1 : ZeroPadPoints) ; %0 pad past pr2*nu2 
for PulseCounter= (dly+1) ;pr2 : (pr2*nu2+ (dly+1) ) ; %pr2*nu2=Nmbr 
for CutUpPulse=0 ;pr2; %...pts in Del. Pulse Train 

if CutUpPulse<=wi2 

DelayedPulseTrain (PulseCounter+CutUpPulse) =1.0; 

end 

end 

end 



if option==2 

%%%%%%% Create Channel 2 Delayed Pulse Train %%%%%%%%%%%%%% 
DelayedPulseTrain2=zeros (1 : ZeroPadPoints) ; %0 padpast pr2*nu2 
for PulseCounter= (dlyCh2+l) :prCh2 : (prCh2*nuCh2+ (dlyCh2+l) ) ; 
for CutUpPulse=0 ; prCh2 ; 
if CutUpPulse<=wiCh2 

DelayedPulseTrain2 (PulseCounter+CutUpPulse) =1 . 0; 
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end 



end 

end 

end 



clc 

echo on 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Signal to Noise Ratio Calculation 

% Define Rayleigh noise power to equal 2. 

% Rician pdf created for sinusoid + Gaussian noise. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
echo off 

SNR=input(' Desired SNR ; ')/ 

NoisePwr=2/ %rayleigh noise pwr =2 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Sum Channel 1 & 2 signals plus their initial delay. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SigEnergy=0; %initialize Channel 1 energy to 0 

for index=l :NumberPointsDel / %sum over Ch. 1 delayed signal 
SigEnergy=SigEnergy+ ( (DelayedPulseTrain (index) ) . '' 2 ) ; 

end 

if option==2 

SigEnergy2=0; %initialize Channel 2 energy to 0 

for index=l :NumberPointsDel2/ %sum over Ch. 2 delayed signal 
SigEnergy2=SigEnergy2+ ( (DelayedPulseTrain2 (index) ) . ''2) ; 

end 

end 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Calc. Ch. 1 signal eunplitude to meet required input SNR. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

SNR10=SNR/10; 

SNRLinear=10''SNR10; %SNR in linear units 

SigAmplitude=sqrt (NoisePwr*SNRLinear*NumberPointsDel/SigEner 

gy) ; 



if option“2 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Calc. Ch.2 signal amplitude to meet required input SNR. 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

SNR10=SNR/10; 

SNRLinear=10^SNR10; %SNR in linear units 

SigAmplitude2=sqrt (NoisePwr*SNRLinear*NumberPointsDel2/SigEn 
ergy2) / 
end 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Create Rayleigh Noise matrix for Channel 1. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
rand <' normal ') / %Gaussian mean=0, variance=l 

for index=l :NumberPointsDel, 

GaussNoisel (index) =rand/ %N (0, 1) 

GaussNoise2 (index) =rand/ %N ( 0 , 1 ) 

end 

RayleighNoise=sqrt ( (GaussNoisel) . '' 2 + (GaussNoise2) . ''2) ; 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Create Rayleigh Noise matrix for Channel 2 . 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
i f opt i on==2 

rand (' normal' ) ; %Gaussian mean=0, variance=l 

for index=l : NumberPointsDel2 , 

GaussNoise3 (index) =rand; %N (0, 1) 

GaussNoise4 (index) =rand; %N (0, 1) 

end 

RayleighNoise2=sqrt ( (GaussNoise3) . '' 2 + (GaussNoise4) . ''2) / 
end 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Make Rician amplitude of Channel 1 delayed pulse train 
% and remove the dc component of entire signal. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SE=0; 

for index=l :NumberPointsDel, %1 --> delay + delayed signal 
if DelayedPulseTrain (index) ==1 %if one, make Rician 



DelayedPulseTrain (index) =DelayedPulseTrain (index) . *rician (Si 
gAmplitude) ; 

SE= (DelayedPulseTrain (index) ) . ''2+SE; 
end %end Rician modification loop 
end 
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SP=SE/NumberPointsDel ; 



%%%%%%% Add noise to Channel 1 Delayed Pulse Train %%%%%%% 
for index=l :NumberPointsDel, %1 --> delay + delayed signal 
if DelayedPulseTrain ( index) ==0 %no Rayleigh noise to signal 

DelayedPulseTrain (index) =DelayedPulseTrain (index) +RayleighNo 
ise (index) ; 

end 

end 

%%%%% Remove dc component of Channel 1 Pulse Train %%%% 
Chldc=mean (DelayedPulseTrain (1 :NumberPointsDel) ) / 
DelayedPulseTrain (1 ;NuiaberPointsDel) =DelayedPulseTrain (1 :Num 
berPointsDel) -Chide; 



if option==2 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Make Rician amplitude of Channel 2 delayed pulse train 
% and remove dc component of entire signal. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SE=0; 

for index=l : NumberPointsDel2, %1 --> delay + delayed signal 
if DelayedPulseTrain2 (index) ==1 %if one, make Rician 

DelayedPulseTrain2 (index) =DelayedPulseTrain2 (index) . *rician ( 
SigAmplitude2) / 

SE= (DelayedPulseTrain2 (index) ) . ^2+SE; 
end %end Rician modification loop 
end 

SP=SE/NumberPointsDel2 ; 



%%%%%% Add noise to Channel 2 Delayed Pulse Train %%%%%%%% 
for index=l ;NumberPointsDel2, %1 --> delay + delayed signal 
if DelayedPulseTrain2 (index) ==0 %no Rayleigh noise to sig 

DelayedPulseTrain2 (index) =DelayedPulseTrain2 (index) +Rayleigh 
Noise2 (index) ; 

end 

end 

%%%%% Remove dc component of Channel 2 Pulse Train %%%% 
Ch2dc=mean (DelayedPulseTrain2 (1 ;NumberPointsDel2) ) / 
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DelayedPulseTrain2 (1 :NumberPointsDel2) =DelayedPulseTrain2 (1 : 

NuiaberPointsDel2 ) -Ch2dc; 

end 

subplot <221) %two rows, two columns 

time= (0 : 1 : ZeroPadPoints-1) ; %Start time axis at 0 

if option==l 

%%%%%%%%%%%%% Plot Reference Pulse Train %%%%%%%%%%%%%%%%%% 

plot (time, ReferencePulse) ; 

title (' Reference Pulse Sequence')/ 

xlabel ( ' time' ) / 

grid 

end 



%%%% Plot Channel 1 pulse train + Rayleigh/Rician noise%%%% 

topDel=l . l*max (DelayedPulseTrain) ; 

bottomDel=l . l*min (DelayedPulseTrain) ; 

axis([0 NumberPointsDel bottomDel topDel]) 

plot (time, DelayedPulseTrain) ; 

title (' Channel 1 Pulse Sequence + Rician/Rayleigh Noise')/ 

xlabel ( ' time' ) / 

grid 



if option==2 

%%%% Plot Channel 2 pulse train + Rayleigh/Rician noise %%% 
topDel2=l . 1 *max (DelayedPulseTrain2) / %10 percent headroom 
bottomDel2=l . l*min (DelayedPulseTrain2) / 
axis([0 NumberPointsDel2 bottomDel2 topDel2]) 
plot (time, DelayedPulseTrain2) / 

title (' Channel 2 Pulse Sequence + Rician/Rayleigh Noise')/ 

X 1 abe 1 ( ' t ime ' ) / 

grid 

end 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Correlation 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

if option==l 

fR=fft (ReferencePulse) / 

fD=fft (DelayedPulseTrain) / 

end 
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if option==2 

fR=fft (DelayedPulseTrain) ; 
fD=fft (DelayedPulseTrain2) ; 
end 

2 =fD.*conj (fR) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Calculate normalizing coefficient 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
Norml=0; 

Norm2=0; 

for k=l : 1 : length (fD) 

Norml= (abs (fD (k) ) . ''2) +Norml; 

Norm2= (abs (fR(k) ) . ''2) +Norm2; 

end 

Norml=sqrt (Norml) / 

Norm2=sqrt (Norm2) / 

NormCoef f =Norml *Norm2 / 
z=z . /NormCoef f; 

z=z*length (fD) ; 

ifftz=ifft (z) ; 
magz=abs (ifftz) / 



%%%%% Flip vector for appearance %%%%%%%%%%%%%%% 
Topmagz= (length (magz) /2+1) / 

Bottommagz=length (magz) /2/ 

magzLen=length (magz) ; %measure length of magz 
Transmagz= [magz (Topmagz tmagzLen) ,magz (1 :Bottommagz) ] 

%%%%%% Make -time to +time axis %%%%%%%%%%%%%%%% 
minustime=- (length (Transmagz) /2) ; 
posittime= (length (Transmagz) /2) -1 ; 
timel= (minustime ; 1 ;posittime) ; 

topmagz=*l . l*max (magz) ; %add 10 percent head room 
axis ( [minustime posittime 0 topmagz]) %scale axis 

subplot (212) 

plot (timel , Transmagz) ; 

title('FFT Correlation : Channel 1 vs. Channel 2')/ 
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xlabel ( ' time lag')/ 
grid 

gtext('+7dB SNR') 

lj4print %Spanagel Room 427 Laserjet 

axis([l 2 3 4])/ %reset axis scaling 
axis/ 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Program Name : CorrPfaLoop.m 19 July 91 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Purpose : To compute the Pfa for Signal / Noise only for 
% then correlation function. This program creates zero mean 
% RAYLEIGH noise. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear 

clg 

subplot (221) 
rand ( ' normal ' ) 

PFA=1.0/ %will be divided by ten to start 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
for PFAC0UNTER=1 : 1 ; 4 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

PFA=PFA/10/ %Pfa = 0.1 0.01, 0.001, 0.001 

%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

for SNR=0:1:18 %Outer loops through all SNR's. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

Iterations=30/ %30 data points per SNR point 

%%%%%%%%%%%%%%%%%%%%%%%%%%% 

for loop=l ; 1 ; Iterations %inner loop, creates groups of +- a 

%%%%%%%%%%%%%%%%%%%%%%%%%%% 

loop/ 

AboveThreshold (loop) =0/ 

BelowThreshold (loop) =0/ 
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Figure 24 . 



FreqCorr7.in MATLAB flowchart. 
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Threshold (loop) =0 ; 

MaxRxy (loop) =0/ 

RxyThreelag (loop) =0; 

%%%%%%%%%%%% Set Parameters %%%%%%%%%%%%%%%%%%%%%%%% 

N=100; %100 noise points, therefore 200 correlation points 

a=10; %+- 10 % on either side of 0 time lag 

Pfa=PFA; %Reestablish because Pfa destroyed each loop 
PfaDef ined=Pf a; 

%%%%%%%%%%%%% Ma)ce N points Noise %%%%%%%%%%%%%%%%%% 
snrnoiseloop %call signal plus noise generation at a SNR 



%%%%% Correlation of zero mean noise only for Pfa calcul . %% 
MeanRayl=mean (RayleighNoise) ; %noise created in snrnoise.m 
MeanRay2=mean (RayleighNoise2) ; 
RayleighNoise=RayleighNoise-MeanRayl ; 
RayleighNoise2=RayleighNoise2-MeanRay2; 

Rxy=xcorr (RayleighNoise, RayleighNoise2) ; %xcorr signal+noise 
later 

%%%%%%%%%%%%%%%%%% Compute t~ %%%%%%%%%%%%%%%%%%%%%%% 
t=0; 

for i=(N-a) :1: (N+a) 
t=t+Rxy (i) ; 

end 

t=t/ (2*a+l) ; 

%%%%%%%%%%%%%%% Compute variance t %%%%%%%%%%%%%%% 

Vart=0; 

for i=(N-a) ;1; (N+a) 

Vart=Vart+ (Rxy (i) -t) .'' 2 ; 

end 

Vart=Vart/ (2*a+l) ; 

%%%%%%%%%%%%%%%%% Define Pfa %%%%%%%%%%%%%%%%%%%% 

Pfa=Pfa*2; %multiply by 2 

Pfa=l-Pfa/ %subtract one 

Threshold (loop) =inverf (Pfa) / 

Threshold (loop) =Threshold (loop) *sqrt (2) *sqrt (Vart) / %mult by 
the sqrt (2) 

%%%%% Xcorr of signal + zero mean noise %%%%%%%%%% 

Rxy=xcorr (DelayedPulseTrain, DelayedPulseTrain2) ; 
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%%%%%%% Compare +a to -a to Threshold %%%%%%%%%%%% 
for i=(N-a) :1: (N+a) 

if Rxy (i) >=Threshold (loop) 

AboveThreshold (loop) =AboveThreshold (loop) +1 ; 
else 

BelowThreshold (loop) =BelowThreshold (loop) +1 ; 

end 

if Rxy (i) >MaxRxy (loop) 

MaxRxy (loop) =Rxy (i) / 

SaveThreshold (loop) =Threshold (loop) ; 

%SaveVart (loop) =sqrt (Vart) ; 

end 

end 

RxyThreelag (loop) =Rxy (97) ; 

end %end # Iterations loop 

%%%%%%% Count correct threshold crossings %%%%%%%% 

CorrectTDOA=0 / 

for count=l : 1 : Iterations 

if RxyThreelag (count) ==MaxRxy (count) 

Correct TDOA=Correct TDOA+1 / 

end 

end 

PercentCorrectTDOA= (CorrectTDOA*100) / Iterations; 

diary flag4 

PercentCorrectTDOA 

BelowThreshold 

AboveThreshold 

PfaDefined 

SNR 

SaveThreshold %compare threshold to zeroth 

lag below 

RxyThreelag %Channel 2 lags Channel 1 by 

three pulses. 

MaxRxy 
diary off 



end 


%end 


# test points loop 


end 


%end 


Pfa loop 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Program Name : snrnoiseloop .m 19 JULY 91 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

dly=0; %Channel 1 delay 

wi2=5; %Channel pulse width 

pr2=10; %Channel 1 pulse period 

nu2=10; %Channel 1 number of pulses 

NumberPointsDel= <pr2*nu2) +dly ; %calc. nmbr pts sig + delay 



dlyCh2=3; %Channel 2 delay 

wiCh2=5; %Channel 2 pulse width 

prCh2=10; %Channel 2 pulse period 

nuCh2=10; %Channel 2 number of pulses 

NumberPointsDel2= (prCh2*nuCh2) +dlyCh2; %nmbr pts sig + delay 



ZeroPadPoints=NumberPointsDel2+NumberPointsDel-l; %N3=N2+N1-1 
for m=2:l:19; %2''(19) = 524280 point FFT max 

if 2‘'m >= ZeroPadPoints, 

ZeroPadPoints=2''m; 

break; %when find a 2^m power, exit loop 

end 

end 



%%%% Create Channel 1 Delayed Pulse Train %%%%%%%%%%% 
DelayedPulseTrain=zeros (1 : ZeroPadPoints) / %0 pad past pr2*nu2 
for PulseCounter= <dly+l) ;pr2: (pr2*nu2+ (dly+1) ) / %pr2*nu2=Nmbr 
for CutUpPulse=0 :pr2; %...pts in Del. Pulse Train 

if CutUpPulse<=wi2 

DelayedPulseTrain (PulseCounter+CutUpPulse) =1.0; 

end 

end 

end 

%%%%% Create Channel 2 Delayed Pulse Train %%%%%%%%%% 
DelayedPulseTrain2=zeros (1 : ZeroPadPoints) ; %0 padpast pr2*nu2 
for PulseCounter= (dlyCh2+l) :prCh2: (prCh2*nuCh2+ (dlyCh2+l) ) ; 
for CutUpPulse=0 :prCh2; 
if CutUpPulse<=wiCh2 

DelayedPulseTrain2 (PulseCounter+CutUpPulse) =1.0; 

end 

end 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Signal to Noise Ratio Calculation 

% Define Rayleigh noise power to equal 2. 
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% Rician pdf created for sinusoid + Gaussian noise. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
echo off 

NoisePwr=2; %rayleigh noise pwr =2 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Sum Channel 1 & 2 signals plus their initial delay. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SigEnergy=0/ %initialize Channel 1 energy to 0 

for index=l :NumberPointsDel/ %sum over Ch. 1 delayed signal 
SigEnergy=SigEnergy+ ( (DelayedPulseTrain (index) ) . ^ 2 ) ; 

end 

SigEnergy2=0 ; %initialize Channel 2 energy to 0 

for index=l :NumberPointsDel2; %sum over Ch. 2 delayed signal 
SigEnergy2=SigEnergy2+ ( (DelayedPulseTrain2 (index) ) . ^2) ; 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Calc. Ch. 1 signal amplitude to meet required input SNR. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

SNR10=SNR/10; 

SNRLinear=10''SNR10/ %SNR in linear units 

SigAmplitude=sqrt (NoisePwr*SNRLinear*NumberPointsDel/ SigEner 

gy) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Calc. Ch.2 signal amplitude to meet required input SNR. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

SNR10=SNR/10; 

SNRLinear=10^SNR10 ; %SNR in linear units 

SigAmplitude2=sqrt (NoisePwr*SNRLinear*NumberPointsDel2/ SigEn 
ergy2) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Create Rayleigh Noise matrix for Channel 1. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
for index=l :NumberPointsDel , 

GaussNoisel (index) =rand; %N (0, 1) 

GaussNoise2 (index) =rand; %N (0, 1) 

end 

RayleighNoise=sqrt ( (GaussNoisel) . ''2+ (GaussNoise2) . '' 2 ) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Create Rayleigh Noise matrix for Channel 2. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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%N (0, 1) 
%N(0, 1) 



for index=l : NumberPointsDel2 , 

GaussNoiseS (index) =rand; 

GaussNoise4 (index) =rand; 

end 

RayleighNoise2=sqrt ( (GaussNoise3) . ''2+ (GaussNoise4) .^2); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Make Rician amplitude of Channel 1 delayed pulse train 
% and remove the dc component of entire signal. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SE=0/ 

for index=l : NumberPointsDel, %1 --> delay + delayed signal 
if DelayedPulseTrain (index) ==1 %if one, make Rician 

DelayedPulseTrain (index) =DelayedPulseTrain (index) . *rician (Si 
gAmplitude) ; 

SE= (DelayedPulseTrain (index) ) . '^2+SE/ 
end %end Rician modification loop 

end 

SP=SE/NumberPointsDel ; 

%%%%%%% Add noise to Channel 1 Delayed Pulse Train %%%%%%% 
for index=l : NumberPointsDel, %1 --> delay + delayed signal 
if DelayedPulseTrain (index) ==0 %no Rayleigh noise to signal 

DelayedPulseTrain (index) =DelayedPulseTrain (index) +RayleighNo 
ise (index) / 
end 
end 

%%%%% Remove dc component of Channel 1 Pulse Train %%%% 
Chldc=mean (DelayedPulseTrain (1 : NumberPointsDel) ) ; 
DelayedPulseTrain (1 : NumberPointsDel) =DelayedPulseTrain (1 ;Num 
berPointsDel) -Chide; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Make Rician amplitude of Channel 2 delayed pulse train 
% and remove dc component of entire signal. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SE=0 ; 

for index=l :N\amberPointsDel2 , %1 --> delay + delayed signal 

if DelayedPulseTrain2 (index) ==1 %if one, make Rician 

DelayedPulseTrain2 (index) =DelayedPulseTrain2 (index) . *rician ( 
SigAmplitude2) / 

SE= (DelayedPulseTrain2 (index) ) . ''2+SE; 
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end %end Rician modification loop 
end 

SP=SE/NumberPointsDel2/ 

%%%%%% Add nc 3e to Channel 2 Delayed Pulse Train %%%%%%%% 
for index=l :NumberPointsDel2, %1 --> delay + delayed signal 
if DelayedPulseTrain2 (index) ==0 %no Rayleigh noise to sig 

DelayedPulseTrain2 (index) =DelayedPulseTrain2 (index) tRayleigh 
Noise2 (index) ; 

end 

end 

%%%%% Remove dc component of Channel 2 Pulse Train %%%% 
Ch2dc=mean (DelayedPulseTrain2 ( 1 : NumberPointsDel2) ) / 
DelayedPulseTrain2 ( 1 : NumberPoint sDel2 ) =DelayedPulseTrain2 ( 1 : 
NumberPointsDel2) -Ch2dc; 

DelayedPulseTrain=DelayedPulseTrain (1 : 100) ; 
DelayedPulseTrain2=DelayedPulseTrain2 (1 : 100) ; 

RayleighNoise=RayleighNoise (1 : 100) / 
RayleighNoise2=RayleighNoise2 (1:100) / 
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Time Domain Correlation 



in Noise Flowchart 




Figur* 25. CorrPfaLoop.m MATLAB flowchart. 
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APPENDIX C. TDOA CONSTANT THKESHOLD SIMULATION 



A. INTRODUCTION 

Four graphs of the output of a time domain based 
correlation detector using a constant threshold are given. 
Each figure has four simulations using a SNR of IdB. The 
is varied from 0.01 to 0.00001 from Figure 26 to Figure 29. 
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lime lag time lag 
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time lag time lag 
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APPENDIX D. FFT OUTPUT NOISE TO SIGNAL RATIO 



A. INTRODUCTION 

The noise to signal ratio for thirteen passes through a 
FFT butterfly are calculated below. The calculation assumes 
the word length of the arithmetic logic unit (ALU) of the 
computer is 16 bits long. The length of the input data points 
are 12 bits long. The data points are right justified when 
placed into the ALU. For a 16 bit ALU and a 12 bit data 
input, scaling occurs at the seventh pass and every other pass 
after that. 

Pass 1 

The noise power at the output of the butterfly after the 
first pass N^^(l) is a function of 

iVjd) = 2 nU 0 ) +205(0) + S2ot(0) + o|(0) 

where = signal power input Co butterfly 

o\ - truncation noise power (74) 

o| = scaling noise power 
oi = sine/ cosine table noise power 
N^{ 0 ) = noise power input to FFT . 

During pass one there is no scaling, no truncation, and the 
trig table is not used. These error terms are zero. Equation 
74 reduces to 
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( 75 ) 



/vi ( 1 ) = 2nI 

where 

= noise input to butterfly 



Pass 2 

The noise power at the output of the butterfly after the 
second pass N^^(2) is a function of 



nI(2) = 2Na(1) + 20s( 1) + 2S2a^(l) +a|(l) . (76) 

During pass two scaling and truncation are not performed, and 
the trig table is not used. These noise terms are zero. 
Substituting Equation 75 into 76 



nI ( 2 ) = 2 {2N^) 

= . 



(77) 



Pass 3 

The noise power at the output of the butterfly after the 
third pass N^^(3) is a function of 

ivi(3) =2Ni{2) + 2o|(2) * 2^S^oi(2) + 0 ^( 2 ) . (78) 

Scaling noise is zero because it is not performed during pass 
three. Substituting Equation 77 into 78 
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Wj(3) = 2[2^nI) + 2^S^qI + o\ 
= S ^ nI + 225^0^ + Q% . 



( 79 ) 



Pass 4 

The noise power at the output of the butterfly after the 
fourth pass N,^^(4) is a function of 

iVi(4) = 2Wj(3) + 2o|(3) + 2^s2o^(3) +0^(3) . (80) 

Scaling noise is zero because it is not performed during pass 
four. Substituting Equation 79 into 80 

nIU) = 2Vg + 2^S^ol + 30r . (81) 



Pass 5 

The noise power at the output of the butterfly after the 
fifth pass N;^^(5) is a function of 

Wj(5) =2iVi(4) + 2o|(4) + 2^S2ot(4) +0^(4) . (82) 

Scaling noise is zero because it is not performed during pass 
five. Substituting Equation 81 into 82 

Ni(5) = 2^nI + (2® + 2^)S^<jI + 70r . (®3) 



Pass 6 

The noise power at the output of the butterfly after the 
sixth pass Nj^^(6) is a function of 
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Wi(6) = 2Na{5) + 2o|(5) ^ 2-S^ai{5) +0^(5) . (8^) 

Scaling noise is zero because it is not performed during pass 
six. Substituting Equation 83 into 84 

iv|(6) = + 2'’S^al + 150 r . (85) 



Pass 7 

The noise power at the output of the butterfly after the 
seventh pass N;^^(7) is a function of 

NaH) = 2 w |( 6 ) + 2 o |( 6 ) + 2^S^ol{6) + 0 ^( 6 ) . ( 86 ) 

Scaling is performed for the first time during pass seven. 
Substituting Equation 85 into 86 

N^il) = 2''Ng + (2® + 2^)S^ol + 310r + 2al . (87) 



Pass 8 

The noise power at the output of the butterfly after the 
eighth pass N,^^(8) is a function of 

Wj(8) =2Nj(7) + 2o|(7) 2^S^ain) +0^(7) . (88) 

Scaling is not performed during this pass. Only the scaling 
noise from the previous pass is added. Substituting Equation 
87 into 88 

Ni(8) = 2®Ng + (2® + 2^)S^al + 630r + 4o| . (89) 
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Pass 9 



The noise power at the output of the butterfly after the 
ninth pass N;^^(9) is a function of 

nI(9) =2Wj(8) + 2o|(8) + 2®S2 o^(8) +0^(8) . (90) 

Substituting Equation 89 into 90 

iv|(9) = 2^Ng + (2^° + 2® + 2^)S^ol + 1270^ + lOo^ . (91) 



Pass 10 

The noise power at the output of the butterfly after the 
tenth pass N,^^(10) is a function of 

iv|(10) =2wi(9) + 2o|(9) + 2®52ot(9) +0^(9) . (92) 

Scaling is not performed this pass. Only the scaling noise 
from the previous pass is added. Substituting Equation 91 into 
92 

wjdO) = 2^°Ng + 2 ^ 25 ^ 0 ^ + 255o| + 2Qcl . (93) 

Equation 93 describes the noise power at the output of a 
butterfly after the tenth data pass. The noise to signal ratio 
after the tenth pass through the butterfly is expressed as the 
ratio 
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( 94 ) 



nI ^ ^ 2^°Ng * 2^^S^ol + 2550^ + 20ol 

_d(10 = H 1 i £ 

52 21 °S 2 



Dividing through by the denominator (signal power) each 
independent and uncorrelated noise term can be identified 






( 10 ) 



where 



Na 



(10) 



52 

4ol 

„2 

4 52 



51 S‘ 




4ot + 4 



_2 

— 

4 52 



51 52 



= output Noise to Signal ratio 

= input Noise to Signal ratio 
= sine/ cosine table noise 
= truncation Noise to Signal ratio 

= scaling Noise to Signal ratio . 
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The noise to signal ratio calculation for 13 passes ( 2 ^^ 
= 8192 points) through the butterfly algorithm continues 



Pass 11 

The noise power at the output of the butterfly after the 
11*^’' pass N*2(11) is a function of 

nUiI) =2Nj(10) + 2c|(10) + 2^®S2o»(10) + 0r(lO) • (^6) 

Substituting Equation 93 into 96 

Nl(ll) = 2“W| + (2^2 ^ * 5110^ + 42c| . (97) 
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Pass 12 



The noise power at the output of the butterfly after the 
12 '^*' pass N;^^( 12 ) is a function of 

Wi( 12 ) = 2 Ni(ll) + 2 (Js(ll) + 2 “sVw(ll) + a'r(ll) . (98) 

Scaling is not performed this pass. Only scaling noise from 
the previous pass is added. Substituting Equation 97 into 98 

nU12) = + (2^* + 2^^)S^al + 1023aT + 84<j| . (99) 



Pass 13 

The noise power at the output of the butterfly after the 
IS^** pass N,^^(13) is a function of 

(100) 

Wj(13) = 2N^(12) + 2o|(12) + 2^^S^<J»(12) +aT(12) . 

Substituting Equation 99 into 100 

( 101 ) 

Wi(13) = 2^^nI + (2^5 + 2^3 + 2^^)S^<ji + 2047Ot + 1700s • 

Equation 101 describes the noise power at the output of a 
butterfly after the IS^** data pass. The noise to signal ratio 
after the IS*^** pass is expressed as the ratio 

nI 2 ^^n| + ( 2 ^® + 2 “ + 2 “) + 2047Ot + 170o| (102) 

= 2^552 • 

Dividing through by the denominator, each independent and 
uncorrelated noise term can be identified 
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5.5a^ 



^(13) 



wi 

where (13) 

5 . 5ol 

1 Or 
"4”^ 




^2 _2 
^ 1 Gy ^ 1 

48,2^ 



= output Noise to Signal ratio 

= input Noise to Signal ratio 
= sine/ cosine table noise 
“truncation Noise to Signal ratio 



1 o| 
4W72'^ 



= scaling Noise to Signal ratio . 



(103) 
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